library(tidyverse)
library(magrittr)
library(plotrix)
library(shinystan)
library(tidybayes)
library(ggridges)
library(cmdstanr)
library(loo)
library(furrr)
plan(multisession)
# User-defined functions
source("../function/r/my_functions.R")
# Import data
# We modified a column name because we realized that
# "block" sounds more intuitive than "task" in description of the analysis
df_gamble <- read_csv("../data/preprocessed/gamble.csv") %>%
rename(block = task)
df_solo <- read_csv("../data/preprocessed/solo.csv") %>%
rename(block = task)
df_group <- read_csv("../data/preprocessed/group.csv") %>%
rename(block = task)
df_solo_group <- bind_rows(df_solo, df_group)
df_gabor_solo <- read_csv("../data/preprocessed/solo_gabor.csv") %>%
rename(block = task)
df_gabor_group <- read_csv("../data/preprocessed/group_gabor.csv") %>%
rename(block = task)
df_gabor <- bind_rows(df_gabor_solo, df_gabor_group)
df_questionnaire <- read_csv("../data/preprocessed/questionnaire.csv")
# Import the fitting results
median_powutil <- readRDS("../output/data/gamble/median_powutil.rds")
median_gabor <- readRDS("../output/data/gabor/median_gabor.rds")
# cmdstanr
# install_cmdstan(version = "2.28.0")
Now we have the two perceptual parameters (\(\epsilon\) and \(\lambda\)) for each participant. The expected probability of having made a correct decision (\(d\)) is expressed as follows (Navajas et al., 2017):
\[ \hat{p}(\mathrm{correct}|\bar{\mu}_{30}, \sigma_{30}, d) = \int_{-\infty}^{+\infty}{d\mu_{30}\hat{p}(\mathrm{correct|\mu_{30}, \sigma_{30}})p(\mu_{30}|\bar{\mu}_{30},\sigma_{30},d)}, \]
The first term inside integral indicates the probability of the correct response given \(\mu_{30}\) and \(\sigma_{30}\). The second term indicates the probability density of \(\mu_{30}\) under \(\mathcal{N}(\bar{\mu}_{30}, \sigma_{30})\).
Firstly, let us review the model.
The participant’s estimate of \(\mu_{30}\) converges to a normal distribution below:
n_theta <- 30
min_theta <- -13
max_theta <- 19
epsilon <- 1
lambda <- 0.95
-50:150 %>%
map(~stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
)) %>%
imap_dfr(~list(
p = dnorm(.y - 51, mean = .x$mu, sd = .x$sigma),
x = .y - 51
)) %>%
mutate(positive = (x > 0)) %>%
ggplot(aes(x = x, y = p, fill = positive)) +
geom_ribbon(aes(ymin = 0, ymax = p)) +
geom_line() +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
scale_fill_manual(values = c("white", "gray")) +
labs(x = "\u03bc30", y = "Probability density") +
guides(fill = "none")
Mean, SD, and the cumulative probability of \(\mu_{30}\):
stat_updating_model(
n_theta = 30, min_theta = -13, max_theta = 19,
epsilon = 1, lambda = 0.95, seed = 1
) %>%
as_tibble()
Given this distribution, let us assume that the agent judges the orientation is clockwise. Then, how much is the expected probability of the correct response when \(\mu_{30} = 60\)?
mu_30 <- 60
-50:150 %>%
map(~stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
)) %>%
imap_dfr(~list(
p = dnorm(.y - 51, mean = mu_30, sd = .x$sigma),
x = .y - 51
)) %>%
mutate(positive = (x > 0)) %>%
ggplot(aes(x = x, y = p, fill = positive)) +
geom_ribbon(aes(ymin = 0, ymax = p)) +
geom_line() +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
geom_vline(aes(xintercept = mu_30), linetype = "dashed") +
scale_x_continuous(breaks = c(-50, 0, mu_30, 150)) +
scale_fill_manual(values = c("white", "gray")) +
labs(x = "\u03bc30", y = "Probability density") +
guides(fill = "none")
The shaded area (cumulative probability) is:
stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
) %$%
pnorm(mu_30 / sigma)
#> [1] 0.9786262
Apparently, we now have the expected probability of the correct response when \(\mu_{30} = 60\).
However, we need to multiply it by the conditional probability density of \(\mu_{30}\) given the agent’s decision. For example, the probability density of \(\mu_{30} = 60\) given the response is clockwise is:
# Calculate normalization constant
z_const <- stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
) %$%
calc_z_const(mu_bar = mu, sigma = sigma, judge = 1)
z_const
#> [1] 0.9231526
# Probability density given the response is clockwise (+1)
stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
) %$%
dnorm(mu_30, mean = mu, sd = sigma) / z_const
#> [1] 0.01219272
To calculate the expected probability of having made a correct decision, we need to integrate the probability density over \(\mu_{30}\) from \(-\infty\) to \(\infty\).
As Navajas et al. (2017)wrote, the expected probability of having made a correct decision is:
\[ \hat{p}(\mathrm{correct}|\bar{\mu}_{30}, \sigma_{30}, d) = \frac{1}{Z}\int_{-\infty}^{+\infty}d\mu_{30}\frac{e^{-\frac{(\mu_{30} - \bar{\mu}_{30})^2}{2\sigma^2_{30}}}}{\sqrt{2\pi\sigma^2_{30}}}\Theta(\mu_{30}d)\Phi\biggl(\frac{|\mu_{30}|}{\sigma_{30}}\biggr) \]
In the above example, the expected probability of having made a correct decision is the probability under the following curve:
-50:150 %>%
map(~stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
)) %>%
imap_dfr(~list(
p = p_hat_integrand(
mu_30 = .y - 51, mu_bar = .x$mu, sigma = .x$sigma, judge = 1
),
x = .y - 51
)) %>%
ggplot(aes(x, p)) +
geom_ribbon(aes(ymin = 0, ymax = p), fill = "gray") +
geom_line() +
labs(x = "\u03bc30", y = "Probability density")
result_updating_model <- stat_updating_model(
n_theta = n_theta, min_theta = min_theta, max_theta = max_theta,
epsilon = epsilon, lambda = lambda, seed = 1
)
result_updating_model %$%
integrate(
p_hat_integrand, lower = -Inf, upper = Inf,
mu_bar = mu, sigma = sigma, judge = 1
) %>%
.$value / calc_z_const(
mu_bar = result_updating_model$mu,
sigma = result_updating_model$sigma,
judge = 1
)
#> [1] 0.885319
By contrast, if the agent’s decision is anticlockwise, the probability is:
result_updating_model %$%
integrate(
p_hat_integrand, lower = -Inf, upper = Inf,
mu_bar = mu, sigma = sigma, judge = -1
) %>%
.$value / calc_z_const(
mu_bar = result_updating_model$mu,
sigma = result_updating_model$sigma,
judge = -1
)
#> [1] 0.659448
Following the equations above, we calculated the expected probability of a correct response. we call this probability \(\hat{p}\). Note that \(\hat{p}\) ranges from 0.5 to 1 because \(\hat{p}\) must be larger than the chance level, 0.5.
df_gabor_nested <- df_gabor %>%
arrange(block, id, trial) %>%
select(!order) %>%
group_by(block, id, trial) %>%
nest()
p_hat_gabor <- df_solo_group %>%
filter(order == "post" | is.na(order)) %>%
select(block, id, trial, judge) %>%
mutate(judge = if_else(judge == "clockwise", 1, -1)) %>%
left_join(median_gabor, by = "id") %>%
arrange(block, id, trial) %>%
left_join(df_gabor_nested, by = c("block", "id", "trial")) %>%
group_by(block, id, trial) %>%
nest() %>%
mutate(
stat = map(data, ~stat_updating_model(
theta = pluck(.x$data, 1, "ori"),
epsilon = .x$epsilon, lambda = .x$lambda
)),
mu = pluck(stat, 1, "mu"),
sigma = pluck(stat, 1, "sigma")
) %>%
unnest(cols = c(data)) %>%
group_by(block, id, trial) %>%
nest() %>%
mutate(
p_integral = map(data, ~integrate(
p_hat_integrand, lower = -Inf, upper = Inf,
mu_bar = .x$mu, sigma = .x$sigma, judge = .x$judge
)),
z_const = map_dbl(data, ~calc_z_const(
mu_bar = .x$mu, sigma = .x$sigma, judge = .x$judge
)),
p_raw = pluck(p_integral, 1, "value"),
p_hat = p_raw / z_const
) %>%
select(block:trial, z_const:p_hat)
The figure shows mean \(\hat{p} \pm SE\) as a function of task difficulty and correct/wrong responses. This pattern is similar to the result (Fig. 2b) in Navajas et al. (2017).
p_hat_gabor %>%
left_join(df_solo_group, by = c("block", "id", "trial")) %>%
group_by(result, id, ori_var) %>%
summarise(ind_mean = mean(p_hat)) %>%
group_by(result, ori_var) %>%
summarise(mean = mean(ind_mean), se = std.error(ind_mean)) %>%
mutate(ori_var = fct_rev(factor(ori_var))) %>%
rename(Result = result) %>%
ggplot(aes(
ori_var, mean, color = Result, group = Result,
linetype = Result, shape = Result
)) +
geom_point() +
geom_line() +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
width = 0.1, linetype = "solid", show.legend = FALSE
) +
ylim(0.6, 1) +
scale_color_manual(values = c("red3", "#56B4E9")) +
labs(x = "Variance of orientations", y = "Expected accuracy (p-hat)") +
theme(legend.position = "top")
We have obtained the estimates of the expected accuracy of having made a correct response.
However, it does not mean that participants know the exact probability. Rather, participants probably felt the probability in a biased manner. Hereafter, we call this biased probability as subjective probability of having made a correct response.
We then fitted five models on subjective probability to the data.
We first made a data-list for the estimation.
datalist_confidence <- p_hat_gabor %>%
inner_join(df_solo_group, by = c("block", "id", "trial")) %>%
filter(block == "solo") %$%
list(
N = nrow(.),
N_subj = length(levels(factor(id))),
id = id,
choice = as.numeric(choice == "risky"),
reward = payoff_risky + payoff_jitter,
p_hat = p_hat,
rho = median_powutil$rho,
theta_var = ori_var,
scale = 1000
)
See the Stan code (../function/confidence/confidence_accuracy.stan) for more details of the priors.
model_onepar <- cmdstan_model("../function/stan/confidence/confidence_onepar.stan")
fit_onepar <- model_onepar$sample(
datalist_confidence, seed = 1, refresh = 100,
chains = 4, iter_warmup = 500, iter_sampling = 2000
)
fit_onepar_rstan <- cmdstan2rstan(fit_onepar)
launch_shinystan(fit_onepar_rstan)
model_twopar <- cmdstan_model("../function/stan/confidence/confidence_twopar.stan")
fit_twopar <- model_twopar$sample(
datalist_confidence, seed = 1, refresh = 100,
chains = 4, iter_warmup = 500, iter_sampling = 2000
)
fit_twopar_rstan <- cmdstan2rstan(fit_twopar)
launch_shinystan(fit_twopar_rstan)
model_prospect <- cmdstan_model("../function/stan/confidence/subjective_accuracy/confidence_prospect.stan")
set.seed(1)
fit_prospect <- model_prospect$sample(
datalist_confidence, seed = 1, refresh = 100,
chains = 4, iter_warmup = 500, iter_sampling = 2000,
init = function(){list(
alpha = runif(63, 20, 30),
tau = runif(63, 10000, 15000)
)}
)
fit_prospect_rstan <- cmdstan2rstan(fit_prospect)
launch_shinystan(fit_prospect_rstan)
model_goldstein <- cmdstan_model("../function/stan/confidence/subjective_accuracy/confidence_goldstein.stan")
fit_goldstein <- model_goldstein$sample(
datalist_confidence, seed = 1, refresh = 100,
chains = 4, iter_warmup = 500, iter_sampling = 2000
)
fit_goldstein_rstan <- cmdstan2rstan(fit_goldstein)
launch_shinystan(fit_goldstein_rstan)
It is difficult to estimate the parameters, we first estimated the parameters for each participant using the L-BFGS method.
model_wugonzalez_solo <- cmdstan_model("../function/stan/confidence/subjective_accuracy/confidence_wugonzalez_solo.stan")
set.seed(1)
fit_wugonzalez_solo <- p_hat_gabor %>%
inner_join(df_solo_group, by = c("block", "id", "trial")) %>%
left_join(median_powutil, by = "id") %>%
filter(block == "solo") %>%
group_by(id) %>%
nest() %>%
mutate(
fit = map(data, ~optimize_confidence_function(
.x, model = model_wugonzalez_solo, n_iter = 100,
refresh = 0, algorithm = "lbfgs", tol_rel_grad = 1e+18, iter = 2000,
)),
estimates = map(fit, ~{.x$mle()}),
lp = map_dbl(fit, ~{.x$lp()})
) %>%
unnest(cols = c(estimates)) %>%
mutate(param = case_when(
row_number() %% 3 == 1 ~ "tau",
row_number() %% 3 == 2 ~ "alpha",
row_number() %% 3 == 0 ~ "beta",
)) %>%
pivot_wider(names_from = param, values_from = estimates) %>%
group_by(id) %>%
mutate(max_lp = max(lp)) %>%
filter(lp == max_lp)
We then performed the hierarchical Bayesian analysis, specifying the initial values.
model_wugonzalez <- cmdstan_model("../function/stan/confidence/subjective_accuracy/confidence_wugonzalez.stan")
fit_wugonzalez <- model_wugonzalez$sample(
datalist_confidence, seed = 1, refresh = 100,
chains = 4, iter_warmup = 1000, iter_sampling = 2000,
init = function() {
fit_wugonzalez_solo %>%
mutate(
alpha = if_else(alpha < 2 & beta < 2, 5, alpha),
beta = if_else(alpha == 5 & beta < 2, 0.1, beta),
alpha = if_else(id == 55, 7, alpha),
beta = if_else(id == 55, 0.5, beta),
alpha = if_else(id %in% c(4, 15, 44), 10, alpha),
beta = if_else(id %in% c(4, 15, 44), 0.05, beta)
) %$%
list(alpha = alpha, beta = beta, tau = tau)
}
)
fit_wugonzalez_rstan <- cmdstan2rstan(fit_wugonzalez)
launch_shinystan(fit_wugonzalez_rstan)
We compared these models using approximate LOO-CV. The Goldstein-Einhorn function outperformed the other models.
loo_confidence <- my_loo(
fit_onepar, fit_twopar, fit_prospect,
fit_goldstein, fit_wugonzalez
)
loo_compare(loo_confidence$loo)
#> elpd_diff se_diff
#> model4 0.0 0.0
#> model2 -8.9 3.1
#> model3 -26.2 5.4
#> model5 -36.1 8.9
#> model1 -248.9 22.9
The result remained unchanged when using WAIC.
waic_confidence <- my_loo(
fit_onepar, fit_twopar, fit_prospect,
fit_goldstein, fit_wugonzalez, method = "waic"
)
loo_compare(waic_confidence$waic)
#> elpd_diff se_diff
#> model4 0.0 0.0
#> model2 -8.9 3.1
#> model3 -26.3 5.4
#> model5 -36.2 8.9
#> model1 -249.0 22.9
The posterior samples cover the actual choice frequency.
fit_goldstein %>%
spread_draws(y_pred[id][ori_var]) %>%
sample_draws(ndraws = 20, seed = 1) %>%
ungroup() %>%
mutate(
n_risky = y_pred,
draw = factor(.draw),
ori_var = case_when(
ori_var == 1 ~ 8,
ori_var == 2 ~ 16,
ori_var == 3 ~ 32,
ori_var == 4 ~ 64
)
) %>%
ggplot(aes(n_risky, ..density..)) +
geom_line(aes(group = draw), stat = "density", color = "gray", size = 0.5) +
geom_line(
data = df_solo_group %>%
filter(block == "solo") %>%
group_by(id, ori_var) %>%
summarise(n_risky = sum(choice == "risky")),
stat = "density", color = "#009E73", size = 2
) +
facet_wrap(vars(ori_var)) +
labs(x = "Choice frequency of the risky option", y = "Density") +
theme_facet +
theme(legend.position = "top")
We first saved the two parameters (\(\alpha\) and \(\beta\)) of the Goldstein-Einhorn function as median_goldstein.
\[ q = \frac{\alpha\hat{p}^\beta}{\alpha\hat{p}^\beta+(1-\hat{p})^\beta}, \]
where \(q\) is the subjective probability of a correct response in a trial.
We also calculated the area surrounded by \(y = x\) and the Goldstein-Einhorn function. This is an index of each participant’s over/underconfidence.
median_goldstein <- fit_goldstein %>%
spread_draws(alpha[id], beta[id]) %>%
median_qi() %>%
select(id, alpha, beta) %>%
group_by(id) %>%
nest() %>%
mutate(
q_area_under_p = map(data, ~integrate(
q_integrand, lower = 0.5, upper = 1,
alpha = .x$alpha, beta = .x$beta
)),
q_area_under_p = pluck(q_area_under_p, 1, "value")
) %>%
unnest(data) %>%
ungroup()
median_goldstein
# This is an example of the area surrounded by q and p_hat
param_2 <- median_goldstein %>%
filter(id == 10)
tibble(p = seq(0.5, 1, length = 1000)) %>%
mutate(
alpha = param_2$alpha,
beta = param_2$beta,
q = goldstein_einhorn(p, alpha, beta),
ymin = if_else(p > q, q, p),
ymax = if_else(p > q, p, q),
fill = p > q
) %>%
ggplot(aes(p, q)) +
geom_ribbon(aes(ymin = ymin, ymax = ymax, fill = fill)) +
geom_line() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
scale_fill_viridis_d(direction = -1) +
labs(x = "Objective accuracy", y = "Subjective accuracy") +
guides(fill = "none")
The relationship between \(\hat{p}\) and \(q\) is shown below.
median_goldstein %>%
group_by(id) %>%
nest() %>%
mutate(
q = map(data, ~goldstein_einhorn(seq(0.5, 1, 0.001), .x$alpha, .x$beta))
) %>%
unnest(cols = c("q", "data")) %>%
mutate(p_hat = seq(0.5, 1, 0.001)) %>%
ggplot(aes(p_hat, q, group = id, color = q_area_under_p)) +
geom_line() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
scale_x_continuous(limits = c(0.5, 1)) +
scale_y_continuous(limits = c(0, 1)) +
scale_color_viridis_c() +
labs(x = "Objective accuracy", y = "Subjective accuracy") +
guides(color = "none")
The figures show posterior distributions of each parameter. mutate_yaxis_num() is defined in ../function/r/my_functions.R.
fit_goldstein %>%
spread_draws(alpha[id]) %>%
left_join(mutate_yaxis_num(median_goldstein, column = alpha), by = "id") %>%
ggplot(aes(alpha, factor(yaxis_num))) +
geom_density_ridges(rel_min_height = 0.01) +
scale_y_discrete(breaks = c(seq(10, 60, 10))) +
labs(x = "\u03b1", y = "Participants")
fit_goldstein %>%
spread_draws(beta[id]) %>%
left_join(mutate_yaxis_num(median_goldstein, column = beta), by = "id") %>%
ggplot(aes(beta, factor(yaxis_num))) +
geom_density_ridges(rel_min_height = 0.01) +
scale_x_continuous(breaks = seq(0, 2.5, 0.5)) +
scale_y_discrete(breaks = c(seq(10, 60, 10))) +
labs(x = "\u03b2", y = "Participants")
We saved medians of individual \(q\) as median_q.
median_q <- fit_goldstein %>%
spread_draws(q[i]) %>%
median_qi()
median_q <- p_hat_gabor %>%
inner_join(df_solo_group, by = c("block", "id", "trial")) %>%
filter(block == "solo") %>%
ungroup() %>%
mutate(q = median_q$q) %>%
select(block, id, trial, q, p_hat)
We performed ordered logistic regression predicting confidence rating (min: 1; max: 6) from subjective accuracy.
model_rating <- cmdstan_model("../function/stan/confidence/rating/ordered_logistic_rating.stan")
datalist_rating <- median_q %>%
left_join(df_solo, by = c("block", "id", "trial")) %>%
group_by(id) %>%
mutate(
scaled_q = c(scale(q)),
scaled_p_hat = c(scale(p_hat))
) %>%
ungroup() %$%
list(
N = nrow(.),
L = 6,
N_subj = length(levels(factor(id))),
id = id,
rating = rating,
scaled_q = scaled_q
)
fit_rating <- model_rating$sample(
datalist_rating, seed = 1, refresh = 100,
chains = 4, iter_warmup = 500, iter_sampling = 2000
)
fit_rating_rstan <- cmdstan2rstan(fit_rating)
launch_shinystan(fit_rating_rstan)
We saved point estimates of the parameters as median_rating.
median_rating <- fit_rating %>%
spread_draws(slope[id], cutoff[id][level]) %>%
median_qi() %>%
ungroup() %>%
select(id, level, slope, cutoff)
The figures show posterior distributions of each parameter. mutate_yaxis_num() is defined in ../function/r/my_functions.R.
fit_rating %>%
spread_draws(slope[id]) %>%
left_join(
mutate_yaxis_num(filter(median_rating, level == 1), column = slope),
by = "id"
) %>%
ggplot(aes(x = slope, y = factor(yaxis_num))) +
geom_density_ridges(rel_min_height = 0.01) +
scale_y_discrete(breaks = seq(10, 60, 10)) +
labs(x = "Slope (b)", y = "Participants")
fit_rating %>%
spread_draws(cutoff[id][level]) %>%
ggplot(aes(x = cutoff, y = as.factor(level))) +
geom_density_ridges(rel_min_height = 0.01) +
xlim(-6, 6) +
facet_wrap(vars(id), ncol = 7) +
labs(x = "Cutoff points (c)", y = "Level") +
theme_facet
The green triangle indicates the actual frequency, and the gray line indicates the 95% interval of the posterior predictive samples.
fit_rating %>%
spread_draws(y_pred[id][rating]) %>%
median_qi() %>%
ungroup() %>%
rename(n = y_pred) %>%
ggplot(aes(rating, n)) +
geom_pointinterval(aes(ymin = .lower, ymax = .upper), color = "gray", size = 0.1) +
geom_point(
data = df_solo %>%
group_by(id, rating) %>%
count(),
color = "#009E73", shape = "triangle"
) +
scale_x_continuous(breaks = 1:6) +
facet_wrap(vars(id), ncol = 7) +
labs(x = "Confidence rating", y = "Frequency") +
theme_facet +
theme(panel.grid.major.x = element_blank())